C     MUNILIST   
C     Copyright (C) 1997  Filip Hroch, Masaryk University, Brno, CZ
C                    
C     This program is free software; you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation; either version 2 of the License, or
C     (at your option) any later version.
C     
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License  
C     along with this program; if not, write to the Free Software
C     Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
C
C    E-mail: hroch@physics.muni.cz
C    Web:    http://www.sci.muni.cz/~hroch/index.html
C    Postal address: F.Hroch
C                    Department of Theoretical Physics and Astrophysics
C                    Masaryk University
C                    Kotlarska 2
C                    611 37 Brno
C                    Czech Republic
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      FUNCTION
     *MUNILIST(NFILE,COLUMN,NSTARS,STARS,OUTLST,MIN0,PER,ALPHA,DELTA,
     *         NSUFIX)

      CHARACTER*80 NFILE,FILE,NSUFIX,SUFIX
      CHARACTER*1 FILTER
      INTEGER MUNILIST,I,J,COLUMN,NSTARS,STARS(*),OUTLST
      DOUBLE PRECISION Y,M,D,H,MIN,S,E,DATJD,JD,MIN0,PER,LS,HCOR
      DOUBLE PRECISION ALPHA,DELTA,HELCOR
      REAL X(200),DX(200)

      MUNILIST = -1
      CALL CTOF(NFILE,FILE)
      CALL CTOF(NSUFIX,SUFIX)

C Input FITS file
      OPEN(12,FILE=FILE,STATUS='OLD',ACCESS='DIRECT',RECL=80,IOSTAT=I)
      IF( I.NE.0 ) THEN
        PRINT *,'Error open file ',FILE
        RETURN
      ENDIF
C read date, time and exposure
      CALL RFITS(12,Y,M,D,H,MIN,S,E,FILTER)
C Evaluate Julian date at half of exposure 
      JD = DATJD(Y,M,D + (H + (MIN + S/60.0)/60.0)/24.0)
      JD = JD + E/2.0/86400.0
      CLOSE(12)
C heliocentric correction
      IF( ALPHA.GE.0.0 )THEN
        LS = 279.465 + 0.985647*(JD - DATJD(Y,1d0,1d0))
        HCOR = HELCOR(ALPHA/57.29578,DELTA/57.29578,LS/57.29578)
      ELSE
        HCOR = 0.0
      ENDIF

c What about long filenames and extensions?
c miroslav.broz@usa.net, Dec 10th 1998

      I = INDEX(FILE,'fts') + INDEX(FILE,'FTS') +
     &  INDEX(FILE,'fits') + INDEX(FILE,'FITS')
      IF( I.EQ.0 ) RETURN
      IF( SUFIX.EQ.' ' )THEN
        FILE(I:I+2) = 'MAT'
      ELSE
        FILE(I:) = SUFIX
      ENDIF
      FILE(I+3:I+3) = ' '
C Open magnitude file
      OPEN(13,FILE=FILE,STATUS='OLD',IOSTAT=J)
      IF( J.NE.0 ) RETURN
c      IF( J.NE.0 )THEN
c        FILE(I:I+2) = 'mat'
c        OPEN(13,FILE=FILE,STATUS='OLD',IOSTAT=J)
c        IF( J.NE.0 ) RETURN
c      ENDIF

C Read magnitude data
      CALL RDAT(13,COLUMN,NSTARS,STARS,X,DX)
      CLOSE(13)

C Output
      DO I = 1,NSTARS
        IF( ABS(DX(I)).LT.1E-3 )THEN
          X(I) = 99.999
          DX(I) = 9.999
        ENDIF
      ENDDO         
      
      IF( OUTLST.EQ.1 ) THEN
        WRITE(*,700) INT(Y),INT(M),INT(D),INT(H),INT(MIN),INT(S),
     .               INT(E),FILTER,JD-2.4E6, (X(I),DX(I),I = 1,NSTARS)
      ENDIF
700   FORMAT (1X,I4.4,TR1,5(I2.2,TR1),I3.1,TR1,A1,TR2,F10.4,1X,1000F7.3)
      
      IF( OUTLST.EQ.0 ) THEN
        DO I = 2,NSTARS
          IF( X(1).GT.99.9 .OR. X(I).GT.99.9 )THEN
            X(I) = 99.999
            DX(I) = 9.999
          ELSE
C compute diferences
            X(I) = X(I) - X(1)
C their errors
            DX(I) = SQRT( DX(1)**2 + DX(I)**2 )
          ENDIF
        ENDDO
        IF( PER .LT. 0.0 )THEN
C Output time is in Julian days.
          IF( NSTARS.EQ.1 )THEN
            WRITE(*,710) JD + HCOR - 2.4E6,(X(I),DX(I),I=1,NSTARS)
          ELSE
            WRITE(*,710) JD + HCOR - 2.4E6,(X(I),DX(I),I=2,NSTARS)
          ENDIF
        ELSE
C Phase
          PHASE = (JD + HCOR - 2.4E6 - MIN0)/PER
          PHASE = PHASE - INT(PHASE)
          WRITE(*,710) PHASE,(X(I),DX(I),I=2,NSTARS)
        ENDIF
710     FORMAT(1X,F10.4,1000F7.3)
      ENDIF

      MUNILIST = 0

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE RFITS(NUNIT,Y,M,D,H,MIN,S,E,F)

C Procedura cte FITS soubor a nadefinuje datum, cas a expozici
C a filter
      INTEGER I,NUNIT,J  
      DOUBLE PRECISION Y,M,D,H,MIN,S,E
      CHARACTER*80 RADEK
      CHARACTER*1 F
      
      Y = 0.0
      M = 0.0
      D = 0.0
      H = 0.0
      MIN = 0.0
      S = 0.0
      E = 0.0
      F = ' '
      
      I = 1
30    READ(NUNIT,REC=I,IOSTAT=J) RADEK
        IF( J.NE. 0 ) THEN
          INQUIRE(NUNIT,NAME=RADEK)
          PRINT *,'FITS File=',RADEK(1:60)
          STOP 'Fatal error for input FITS file.'
        ENDIF
        I = I + 1
        IF(RADEK(1:9).EQ.'DATE-OBS=')THEN  
          READ(RADEK(12:20),'(3(E2.0,TR1))') D,M,Y
          Y = Y + 1900.0
        ELSEIF(RADEK(1:9).EQ.'UT      =')THEN
          READ(RADEK(12:20),'(3(E2.0,TR1))') H,MIN,S
        ELSEIF(RADEK(1:9).EQ.'EXPTIME =')THEN
          READ(RADEK(11:30),'(E30.0)') E
        ELSEIF(RADEK(1:9).EQ.'FILTER  =')THEN
          F = RADEK(12:12)
        ELSEIF(RADEK(1:3).EQ.'END')THEN
          GOTO 40
        ENDIF
      GOTO 30
40    CONTINUE

      END
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE RDAT(NUNIT,COLUMN,N,NLIN,X,DX)
      
C  Cteni jasnosti ze souboru. Sam pozna, jestli jde 
C o typ AP nebo jiny a podle toho cte chyby jasnosti.

      INTEGER NUNIT,COLUMN,N,NLIN(*)
      REAL XX,X(*),DX(*)

      CHARACTER*200 LIN
      LOGICAL AP
      INTEGER I,J,IN
      
      DO I = 1,N
        DX(I) = 0.0
      ENDDO

C Rozpoznani typu souboru
      AP = .FALSE.
      READ(NUNIT,'(A)') LIN
      IF( LIN(1:4).EQ.' NL ' )THEN
        READ(NUNIT,*) I
        AP =  I.EQ.2
      ENDIF
      
      REWIND(NUNIT)
      J = 0
           
C Zpracovani radku
50    J = J + 1
        READ(NUNIT,'(A)',END=60) LIN
        IN = 0
        DO I = 1,N
          IF( NLIN(I).EQ.J ) IN = I
        ENDDO
        IF( IN.NE.0 )THEN
          IF( AP )THEN
            READ(LIN,800) (XX,I=1,COLUMN-1),X(IN)
            READ(NUNIT,900) (XX,I=1,COLUMN-1),DX(IN)
            J = J+1
          ELSE
            READ(LIN,800) (XX,I=1,COLUMN-1),X(IN),DX(IN)
          ENDIF
        ENDIF
      GOTO 50
60    CONTINUE         

800   FORMAT(F6.0,14F9.3) 
900   FORMAT(F13.0,2F6.2,F8.3,11F9.3)
                             
      END

